Stability of Leap- Frog Constant-Coefficients Semi-Implicit Schemes for 
the Fully Elastic System of Euler Equations. Flat-Terrain Case. 



P. Benard*, R. Laprise°, J. Vivoda+, P. Smolikova^ 

* Centre National de Recherches Meteorologiques, Meteo-France, Toulouse, France 
° Universite du Quebec a Montreal, Montreal, Canada 
+ Slovak Hydro-Meteorological Institute, Bratislava, Slovakia 
^ Czech Hydro-Meteorological Institute, Prague, Czech Republic 

12 September 2003 

Corresponding address: 

Pierre Benard 

CNRM/GMAP 

42, Avenue G. Coriolis 

F-31057 TOULOUSE CEDEX 

FRANCE 

Telephone: +33 (0)5 61 07 84 63 
Fax: +33 (0)5 61 07 84 53 
e-mail: pierre.benard@meteo.fr 

1 



ABSTRACT 

The stability of Semi-Implicit schemes for the Hydrostatic Primitive Equations system 
has been studied extensively over the past twenty years, since this temporal scheme and 
this system represented a standard for NWP. However, with the increase of computational 
power, the relaxation of the hydrostatic approximation through the use of nonhydrostatic 
fully elastic systems is now emerging for future NWP as an attractive, solution valid 
at any scale. In this context, several models employing the so-called Euler Equations 
together with a constant-coefficients semi-implicit time discretisation have already been 
developed, but no solid justification for the suitability of this algorithmic combination has 
been presented so far, especially from the point of view of robustness. 

The aim of this paper is to investigate the response of this system/scheme in terms 
of stability in presence of explicitly treated residual terms, as it inevitably occurs in 
the reality of NWP. This sudy is restricted to the impact of thermal and baric residual 
terms (metric residual terms linked to the orography are not considered here). It is 
shown that conversely to what occurs with Hydrostatic Primitive Equations, the choice 
of the prognostic variables used to solve the system in time is of primary importance 
for the robustness with Euler Equations. For an optimal choice of prognostic variables, 
unconditionnally stable schemes can be obtained (with respect to the length of the time- 
step), but only for a smaller range of reference states than in the case of Hydrostatic 
Primitive Equations. This study also indicates that: (i) vertical coordinates based on 
geometrical height and on mass behave similarly in terms of stability for the problems 
examined here, and (ii) hybrid coordinates induce an intrinsic instability, the practical 
importance of which is however not completely elucidated in the theoretical context of 
this paper. 
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1 Introduction 



In their most general definition, semi-implicit (SI) schemes consist in an arbitrary separa- 
tion of the evolution terms of any dynamical system between some linear terms, treated 
implicitly, and non-linear residuals (NL residuals hereafter), treated explicitly. In mete- 
orology, depending on the nature of the implicitly-treated terms, three main types of SI 
schemes can be distinguished. The coefficients of these linear terms may be: (i) constant 
both in time and horizontally; (ii) constant in time only; and (iii) non-constant. 

The first approach was initially introduced for meterological applications by Robert 
et al., 1972, and has been extensively used in numerical weather prediction (NWP) since 
the solution of the resulting implicit system requires only basic techniques. However, due 
to large NL residuals, the stability of these schemes is not formally guaranteed, especially 
for long time-steps. 

The second and third approaches require more sophisticated techniques for solving the 
resulting implicit system, but they allow a significant reduction of the magnitude of the 
explicitly treated residuals, and hence, a potentially better stability. 

In the present paper, the terms "constant-coefficients SI schemes" exclusively refer to 
the above first category of SI schemes, and only these schemes are considered in all the 
following, unless expressly mentionned. 

Historically, SI schemes have been first applied in NWP for solving the hydrostatic 
primitive equations (HPE) system, and extensive stability studies have been carried out 
with this system. Simmons et al, 1978 (hereafter SHB78) investigated the practical sta- 
bility of HPEs with the three-time-level (herafter referred to as 3-TL) leap-frog constant- 
coefficients SI scheme in the terrain-following pressure-based a coordinate. To do so, they 

examined the effect of the leading NL residual terms on the stability when the SI reference 

3 



temperature deviates from the actual temperature. They showed that, in the particular 
case where the complete model operator and the linearized SI operator have the same 
eigenfunctions, the stability can be studied analytically. In the more general case, when 
this latter condition is not fulfilled, Cote et al., 1983 (hereafter CBS83) showed that the 
stability can still be assessed, but at the price of a "numerical analysis" which can be 
performed only in the space-discretized context; the stability analysis then becomes an 
eigenvalue problem in a generalized state-vector space where the whole space- and time- 
discretized model acts as a so-called "amplification matrix". The salient result of SHB78 
was that a warm isothermal choice for the SI reference temperature profile resulted in a 
more stable scheme than when using climatological profiles for the SI reference-state, a 
rule which has been widely followed in practical NWP applications. The two methods 
proposed by SHB78 and CBS83 (analysis in simplified cases, and numerical analysis in 
the general case) have been adopted by most of subsequent studies on the SI scheme 
stability. CBS83 showed that, for the finite-element vertical discretization of a 3-TL HPE 
model in a coordinate, the SI reference-state static stability had to be larger than half 
the actual one, in order to achieve stability, which explains and generalizes the previous 
results by SHB78. Simmons and Temperton (1997, hereafter ST97) have extended the 
study of SHB78 to extrapolating two-time-level (hereafter 2-TL) SI schemes, still in the 
HPE system but for the more general hybrid-pressure terrain-following rj coordinate (in 
r] coordinate, the surface pressure in the SI reference-state also has to be considered for 
the stability of the SI scheme). 

However, with the increasing resolutions allowed by faster computers, nonhydrostatic 

NH models are now accessible to NWP, thus avoiding the limitations associated with the 

hydrostatic assumption. However, the aim for NWP models should be to achieve the 

relaxation of the hydrostatic approximation while keeping the same degree of efficiency 
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as the former HPE SI semi-Lagrangian systems. This pleads in favor of attempting to 
extend the SI technique to the new generation of nonhydrostatic semi-Lagrangian models. 
Among the possible NH systems, the fully elastic Euler Equations (EE) system is generally 
advocated for several reasons, including the possibility to build a "universal" atmospheric 
model which dynamical kernel is valid for all scales used by the meteorological community. 

Several models using the EE system with a constant-coefficients SI scheme have been 
recently developed (Tanguay et ah, 1990; Laprise et al, 1995; Semazzi et ah, 1995; Bub- 
nova efc ai., 1995, hereafter BHBG95; Caya and Laprise, 1999), but the stability of such 
a combination has not yet been studied in detail. Tanguay et ah (1990) examined the 
stability of the tangent-linear version of their model around the SI reference state (i.e. in 
the absence of non-linear terms). Not surprisingly, they found that the system is uncondi- 
tionally stable, since this follows from a general property of any purely linear SI physical 
system. However, the experience accumulated with former HPE systems clearly indicates 
that some care must be taken when (explicitly treated) nonlinear terms are present. More- 
over, there is no objective reason to expect an intrinsically more robust behavior of the SI 
EE system compared to the SI HPE system. Especially, when increasing the resolution 
to mesoscales, terms associated with the orographic forcing or with physical processes 
could have an increasingly stringent impact on the stability, due to their increased con- 
tribution in the total evolution. In BHBG95, we reported an unstable behavior for the 
EE system with the SI time discretization, and this problem was solved in the Eulerian 
context through the iteration of a part of the NL residual terms. However the unstable 
behavior reappeared when the resolution was refined or when the timestep was increased, 
as allowed by the implementation of a semi-Lagrangian scheme. This prompted us to 
undertake this study. A general formalism for studying the stability of various time dis- 
cretizations (including the SI scheme) and equation systems (including the EE system) 
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has been presented in Benard (2003), BOS herefater. This general method is applied here 
to the case of the 3-TL SI scheme for the EE system with flat orography, in order to 
show the essential role played by the choice of prognostic variables in the robustness of 
the system. 



2 Framework of analyses 

In this paper (except in Appendix B), the EE are cast in the pure unstretched terrain- 
following coordinate a which can be classically derived from the mass-based hydrostatic- 
pressure coordinate vr (Laprise, 1992) through a = (tt/tTs), where vr^ is the hydrostatic 
surface pressure. The a coordinate examined here is a particular case of the general 
stretched hybrid-pressure terrain-following coordinate r] of BHBG95. However, the use 
of the a coordinate is advantageous for the theoretical analysis since it possesses a much 
simpler vertical metrics. Starting from the general rj formalism, the equations for the cr 
coordinate can be obtained by setting the arbitrary A{ri) and B{ri) functions in equation 
(9) of BHBG95 as follow: 

Mv) = (1) 

B{v) = r/^a (2) 

In all the discussions of this paper, the flow will be assumed adiabatic inviscid and friction- 
less in a non-rotating perfect-gas dry atmosphere with a Cartesian coordinate system. In 
these conditions, the complete set of Euler equations can be easily derived from BHBG95 
equations (l)-(8) with the same standard notations (for non-standard notations see Ap- 
pendix A). The system in a coordinate writes: 
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^ + fir:^ + l*v^ = (3) 

at p TXg oa 

^'■«fl--tl-0 (4) 



dt \ -Kg da 
dT RT 



(VVg + VV)d(T = 



(7) 



where the three-dimensional divergence is given by: 



fl3 = V.V + ^f?^W^-9^^ (8) 



TisRT \da ) ^ " nMT da 
and the geopotential horizontal gradient is: 



V0=V0, + i?/ v(^]da (9) 



P 

We note g = ln(7rs), p is the true pressure, and w is the vertical velocity. 
The theoretical analysis of the stability of non-linear systems such as (jSl)-© is not possible 
in the most general conditions and requires some simplifications in order to becomes 
algebraically tractable. For the analyses presented here, we use the same method and 
notations as in BOS, which is basically a generalisation of the method proposed by SHB78. 
In symbolic notations, the equations of the system to be solved can be written as: 

-Q^ = M{X) (10) 

where X is the state variable, and M. is the full model operator containing only spatial 
dependencies on X . For the definition of the SI scheme, as usually done in NWP, a SI 
reference-state X* is chosen, and the system M. is linearized around X* ^ resulting in a 
linear operator noted £*. The 3-TL SI time discretization writes: 



[Mix'') - c.x^] + r . { ^^^'^^ (11) 



2At 
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where the superscript (-, 0, +) indicates variables at time (t — At), t, and (t + At) 
respectively, where At is the lenth of the timestep. The stability of the model is then 
conditioned by the structure of the NL residual {M. — £*). 

For the purpose of analyses, the flow is assumed to consist of small perturbations 
around a steady basic-state X. Hence A4{X) = (dX/dt) = 0, and the full model evolution 
can be described by C, the linear-tangent operator of around X. 

The 3-TL SI time-discretized equation then becomes: 



which represents the equation to be solved in the simplified framework examined here. 
Note that X X* and hence £ 7^ £* thus giving rise to explicit contributions in (fT2l ) 

The other simplifications adopted for the present analyses (as well as in SHB78), con- 
sist in assuming that both basic- and SI reference-states are resting, isothermal and 
hydrostatically-balanced, with a uniform (plateau) orography. 

As a consequence, in hydrostatic-pressure based coordinate, the X and X* states 
are simply characterized by their uniform temperature and hydrostatic surface-pressure 
fields, i.e. by the pairs of numbers (T, vf^) and (T*, vr*) respectively. Finally, the domain 
is assumed two-dimensional (in a vertical plane) with x as horizontal coordinate. 

The fact to examine the stability of the numerical system in this highly simplified context 
of course leads to an overestimation of the stability compared to what can be expected for 
more complex fiows with a fully non-linear model. However if a scheme is found unstable 
here, it will have little chances to be applicable in practice. 




(12) 
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3 Impact of the Sl-reference pressure 



In this section, we examine the SI scheme for EE when the actual hydrostatic surface- 
pressure vfj deviates from its SI reference counterpart vr*, everything else being equal for 
the two states. Hence, in this section the basic X and reference X* states are characterized 
by (T = T*, Tfs) and (T*, tt*) respectively. The linearized equations are derived from 
(jSI)-© or from Eq. (13)-(21) of BHBG95, with the following set of prognostic variables: 
horizontal divergence D, temperature T, q = In^Tig), and the following two nonhydrostatic 
variables as in BHBG95: 

p — 71 



V 

d 



TT" 



-g- 



Tx* dw 



m*RT* dr] 

where vr*(?7) is the SI reference pressure profile, and m*{vi) 
these variables become: 



V 

d 



p — TT 



an: 



a 



-9 



dw 



RT* da 



(13) 
(14) 

{dn* Idrf). In a coordinate, 

(15) 
(16) 



The C system then writes: 
dD 

'dt ~ ~ 
dd _ _ 
dt ~ ~ 
dT 

'dt ~ ~ 
&P 
'dt 

dq 

di 



-ngV^T + RT* 

g' rvr! 



{g -I)V'V - RT*V'q 



RT* 
RT* 



(D + d) 



d{d + I)V 



"vr/ 




'vr/ 













{D + d) 



-MD 



(17) 

(18) 
(19) 
(20) 
(21) 



where the notations follow BHBG95 (see also Appendix A for the meaning of vertical 

operators X, S, J\f and d in a coordinate). 
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The SI linear model C* can be derived in a similar way except that vr* is used instead 
of Tfs: it is then formally identical to £, except that the factors in brackets in the above 
system become equal to 1. As a consequence, the explicitly-treated NL residual model 
{C — £*) in (|T^ is non-zero, leading to a potential source of instability resulting from the 
departure of vf^ from vr*. 

An alternative formulation can be obtained if the variable V is replaced by a new variable 
V defined by: 

P = ^ (22) 

TT 

The form of the full nonlinear model M. is modified in such a way that all occurences of V 
in the RHS must be replaced by {tx/-k*)V^ and the pressure departure equation becomes: 



dV _ TT* 

dt 71 



dV 



7iV d /vr* 
71* dt \ Ti: 



(23) 



dt 

where the bracketed term is the RHS of the original prognostic V equation. The last term 
of (f23|) is identically zero for the linearized C and C* systems. Moreover the factor {n* /ti) 
writes (7r*/7fs) in the present linearized context in a coordinate. As a consequence, the 
new C system writes: 
dD 

— = -RQW^T + RT*{g -I)W^V - RT*V^q (24) 

dd ~ ~ 

^ d{d + X)V (25) 



dt RT* 
dT RT* 



dt Cy 



{D + d) (26) 



Since this system has no dependency on vr*, it is obvious that the SI reference system C* 

will have exactly the same form as C for this set of variables. In fact, it can be seen that 

using the variables (g, V), and the a coordinate allows to completely eliminate the SI 

reference surface-pressure vr* from the model formulation, in opposition to what occurs 
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when using vr^ and/or V as prognostic variables and/or the general hybrid rj coordinate. 
For the new variable "P, there is no NL residual terms, and hence no potential source of 
instability due to the discrepancy between tXs and tt* for the examined problem. As a 
direct consequence, no stability analysis is necessary here to conclude that the variable V 
is better suited to the design of a SI scheme than V. 

Similar algebraic derivations show that for the particular problem examined here, the 
various possible choices for the prognostic pressure variables fall into two classes: 

(i) variables leading to potentially unstable SI: p, p/poi P ~ ^r, p/tt*, V 

(ii) variables leading to stable SI: ln(p), ln(p/po), p/tt, ln(p/7r), V 

where po is an arbitrary constant. It is worth emphasizing that the above statement 
holds for height-based coordinates as well as for the mass-based coordinate that was used 
here (with of course the restriction due to the fact that the variables involving tt are 
not natural with height-based coordinates). These properties follow immediately from 
the derivation of the corresponding linear system in the same context, for height-based 
coordinates. From now on, the new variable V will be used instead of the original variable 
V used in BHBG95. 

4 Impact of the SI reference temperature 

In this section, we examine the stability of the SI scheme for EE with the prognostic 

variables {D, d, T, V, q), when the basic uniform temperature T deviates from the 

SI reference state temperature T*, everything else being equal for the two states. The 

variable d is still given by (fT6|) . As a consequence, the 3-dimensional divergence (see 

Eq. (20) of BHBG95) writes for the Z system: 
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(29) 



and the direct linearization of the original system yields: 



dD 
'dt 
dd 
di 
dT 

'dt 
&P 
'dt 
dq 
di 




-RQV^T + RT{g - J) V^P - RTV^q 



MD 



(34) 



(31) 



(33) 



(32) 



(30) 



The C* operator is defined in a similar way, simply replacing T by T* in the RHS of the 
above system. 

The method for the stability analysis exactly follows the one proposed in B03, and the 
reader is invited to refer to this paper for more details on the notations and the algebraic 
developments. The above system is first shown to fulfil the four conditions [C1]-[C4] 
required for making possible the space-continuous analyses with the proposed method. 
The number of prognostic variable is P = 4 in the sense of B03, and the space-continuous 
state-vector is A" = (A^i, . . . , .^4) = {D,d,T,V). The linear operator in [CI] involves 
li = d applied to (|3n| and I4 = {d + X) applied to (j33|l as in section 7.1 of B03. The 
condition [C'2] requires T > 0, and the normal modes of the system are then: 



where (/c, i^) G R (note that is a non-dimensional vertical wavenumber). In this partic- 
ular case, the four components (/i, . . . , /4) of the shape function / introduced in B03 are 
identical. The vector function / represents the geometry of any normal mode of the time 
and space-continuous system. The verification of [C3], [C4] proceeds easily, as in B03; for 
[C3], we have: 



Xj{x, a) = Xj exp{tkx) a^*'^^^/^) ^ j e (1, . . . , 4) 



(35) 
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^1 = iv-l/2 
= iu + 1/2, 



(36) 
(37) 



and for [C4]: 



^13 



-k^R 



Jl^^ = k'^RT{iiy + 1/2) , /i*^ = k'^RT*{iv + 1/2) 

-.2 



/^24 
7^31 
7^41 
A''42 



= (z.2 + 1/4) 

RT 
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RT* 

* — * 

/^31 ~ /^32 ~ /^32 



/ill = 1-^(^^ + 1/2) 
-g(zz. + l/2)^, /i^^ 



-^(^^ + 1/2). 



(38) 
(39) 
(40) 
(41) 
(42) 
(43) 



For the stability analysis, the growth of any mode with the shape function / is examined. 
The analysis hence consists in solving (fT2|) assuming X^t=-At) = '^'■fix, a) together with: 



^(t=0) — ^X{t=-At) 



(44) 
(45) 



where the unknowns are the complex polarisation vector X~ and the numerical complex 
growth rate A. As stated in B03, the 3-TL SI scheme is a particular ICI scheme with 
Niter = 1 and //(A) = 2 — 1/A. Hence, in the formalism of B03, the stability problem 
reduces to: 



^ (2A - l)h -h O4 ^ 



V 



.Z = M.Z = 0. 



(46) 



Ml M2 Ms 
-AJ4 O4 h 

where the generalised state-vector Z is defined hy Z = {X~ , 1^ and O4 

are the unit and null 4th-order matrices respectively. In the 3-TL SI framework, the 
sub-matrices Mi, M2 M3 are defined by: 
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(Ml),. 



-At- {fi^j - fi*^) 



(47) 
(48) 
(49) 



where all notations follow BOS. The possible values of A for the normal mode structure 
that we examine, are thus given by the roots of the following polynomial equation in A: 

Det(M) = (50) 

In this simple case, the determinant is easily expanded algebraically, and yields: 



(A- 



+ c2^4^ (A+) ( ^'A + rinA 



AACp - {A+yR 



nH I^AA - (A+)^ 



(51) 



where the time-discretised response factors are defined by: 

A^-l 



A+ 
A 
A 



A2 + 1 



(A+) + a\ 



(A+)- 



a 



1 + a 



-A, 



(52) 
(53) 
(54) 
(55) 



and: 



a 



H 

N 
n 
n 



RT*{CJC,) 
RT*/g 
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(56) 
(57) 
(58) 
(59) 
(60) 
(61) 



This eighth-degree complex polynomial equation in A can be solved numerically: for any 
pair {k, z/), the modulus of the eight roots A give the growth rate of the eight correspond- 
ing eigenmodes (four physical modes and four computational modes, due to the 3-TL 
discretization). If one of the roots has a modulus larger than one, then the corresponding 
mode is unstable. The stability of the scheme for the structure function / correspond- 
ing to a pair {k, u) is then given by the maximum modulus of the eigth corresponding 
eigenvalues: 



r = Max(|Ai|) , « G (1,...,8) (62) 

The criterion for unconditional stability of the scheme with respect to the time-step can 
be found by requiring stability at the large time-steps limit in the above equation. The 
terms containing A_ are then vanishing and the equation for the growth rate becomes: 



(Ah 



(63) 



The four modes represented by (A+)^ = are always neutral. Conversely, in the other 
set of roots the short vertical modes are always unstable. In effect, for these roots, 
substituting the time discretization response factors by their value leads to: 

(1 + a){\^ + if - 2a'X{\ - if [lu + (1/2 - CJC,)] = (64) 



For short modes, u » 1, hence the modulus of the four roots becomes close to 

l + a 



^1(1,2,3,4) ~ TTT.^ (65) 



and obviously this leads to large instabilities. As a consequence, except in the degenerated 
case a = 0, the scheme cannot be unconditionally stable in At, because for large time- 
steps, short enough vertical modes are unstable. 

However, drawing on the results of the previous section, the sensitivity of the stability to 

the choice of the prognostic variables is suspected, and the relevance of the original choice 
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d could be questioned. A close inspection of the algebra in the above analysis indicates 
that the source of the problem lies in the small discrepancy between AA and (A+)^ in 
(|63|) . which in turn is linked to the discrepancy between the D and d factors in the RHS 
of (j^ . This suggests the use of an alternative variable d which would be defined in the 
general hybrid coordinate r] by: 

where m = {dn/drf). In the present linear and hydrostatically-balanced context with cr 
coordinate, d simplifies to: 

d = -g^— 67 

and the linear 3-dimensional divergence writes: 

= ^ + d (68) 

The linear system C becomes: 

BD — — 

— = -RgV^T + RT{g-I)V^V-RTV^q (69) 
dd ~ ~ 

ai = + '™' 

^^(^ + d) ,71, 



dt a 



V 



f = 5fl-|(flH-d) (72) 

Since d and d have the same value in the SI reference state, the general design of the SI 
scheme is unchanged: The linear system C* is formally identical to the previous one, hence 
the modification does not change the SI equation to be solved. The stability analysis for 
this new system can be done exactly in the same way as presented above, and the stability 
equation (f5T| for a given geometry {k, u) becomes: 

+ {^'~^ + ^^^) + k'N'c\K^fll = (74) 
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As previously, the unconditional stability in At can be examined by neglecting all terms 
containing A_, and the numerical growth rate is then given by: 

(A" + 1)\X^ + 1 + 2aA)(A' + 1 - 2^^^) = (75) 

The four modes involved by the first factor are always neutral. Basic algebraic manipula- 
tions show that the modes involved in the second factor are stable for —1 < a < 1, while 
the last factor requires —1/2 < a for stability. Finally, considering the definition of a, 
the scheme is unconditionally stable in At when the following condition is fulfilled: 

— <T <2T* (76) 

The above type of stability analyses can be performed in the same way for any other 
vertical-velocity related prognostic variable. From the set {w , {dw/drj), d, d}, only the 
last variable is then found to allow a non-vanishing range of unconditional stability in 
mass-based coordinates. But conversely to what occurred in the previous section, this 
result is now dependent on the type of vertical coordinate used: for height-based coordi- 
nates, the same kind of analysis (see Appendix B) shows that the two most natural choices 
{w , dw/dz} lead to the same unconditional stability criterion than one obtained here for 
d variable in mass-based coordinate. Hence, from the point of view of the asymptotic 
stability at long time-steps, height-based and mass-based coordinates behave identically 
provided "optimal" variables are chosen. The fact that d in mass-based coordinates and 
(dw/dz) in height-based coordinates behave similarly can be understood intuitively since 
these two variables are in fact two expressions of a same concept in the present simplified 
context (i.e. they both represent the true vertical divergence here). The reason why w 
is stable in height-based coordinates but unstable in mass-based coordinates is more sub- 
tle: as suggested by the above analyses, the robustness of the SI scheme appears to be 
highly correlated to the existence of NL residuals in the elastic term D^; when w is used 

as a prognostic variable in height-based coordinates, can be readily obtained from 
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D + (dw/dz), and has no NL residual, thus leading to a stable scheme in the examined 
context. With mass-based coordinates however, ifw is chosen as a prognostic variable, -D3 
must be evaluated through D — {gp / RT){dw / Dtx) and the existence of a NL residual when 
T ^ T* may lead to instabilities. The fact that the vertical divergence d is completely 
imposed as a prognostic variable if a robust SI scheme is desired for the EE system with 
flat-terrain is thus a speciflcity of mass-based coordinates. 

Another important result from the above analysis is that, even with an optimal choice 
of the prognostic variables, the stability range is dramatically reduced for EE system 
compared to the HPE system, whatever vertical coordinate is used. For HPEs (see e.g. 
SHB78), the criterion (f76|) would write: 

< T < 2r* (77) 

Chosing a very warm T* garantees stability in HPEs while this is no longer the case for 
EEs. Moreover, if the actual temperature were to vary by more than a factor 4 in the 
atmosphere, the SI technique examined here could not offer any stable scheme for the EE 
system. The SI EE system is thus less stable hy nature than the SI HPE system, and 
its applicability for NWP is only made possible thanks to the moderate variability of the 
thermal fleld in the terrestrial atmosphere. 

5 Practical Implications 

For a given mode, the root of maximum modulus in ()64|) and (|75|) gives the asymptotic 

growth rate for large timesteps for the two variables d and d examined above. For the 

variable (i, this asymptotic growth rate is function of a and hence the instability is 

directly linked to the vertical resolution of the model; for the variable d the asymptotic 

growth rate is function of a only. Figure Q] shows the asymptotic growth rate for a 
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vertical spacing of 100 m (i.e. a minimum vertical wavelength of 200 m), for d and d, in 
the conditions of the above analyses. Assuming a typical value ±0.25 for a in realistic 
conditions, the range of growth rate for the variable d is clearly incompatible with a stable 
integration of the SI scheme with very long time-steps. In practical applications however, 
the time-step is bounded, and the growth, if present, may not endanger significantly the 
stability of the scheme. However, direct numerical solution of (f5l|) shows that in fact 
any mode (/c, u) is unstable for any timestep, when a 7^ 0. For instance, in Fig. [T] are 
also plotted the numerical growth rates obtained for a mode which geometry is close to 
the shortest mode for a typical mesoscale future NWP limited-area target configuration 
with 3-TL semi-Lagrangian scheme (Ax = 2000 m, A2; = 100 m, and At = 20 s): the 
system with d variable is significantly unstable as soon as a 7^ 0, even for this modest 
timestep. The d variable is clearly demonstrated as not suitable for use in a SI EE scheme 
in mass-based coordinates. This was one of the reasons why an iterative procedure was 
required and applied in the NH model described in BHBG95. 

An additional remark concerning the stability for finite timesteps is that surprisingly, 
when an optimal variable is chosen in height-based coordinates, the SI scheme which is 
stable for long timesteps becomes slightly unstable for finite timesteps as soon as a 7^ 0, 
as seen in Appendix B. A similar behaviour was found in BOS for the one dimensional 
acoustic system with long time-steps in height-based coordinates. This unstable behaviour 
can be interpreted as originating from the fact that in height-based coordinates, the time- 
continuous normal modes of the C system are not normal modes of £*, as pointed out in 
Appendix B. As a consequence, the time-continuous evolution of any normal mode of C by 
the C* system contains a growing component (i.e. a complex non-real frequency) as soon 
as a 7^ 0. The SI scheme is then not able to ensure a stable evolution for some of these 

components. This type of instability does not occur for mass-based coordinates, as seen 
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in Fig. m From the theoretical point of view, this behaviour is in defavour of height-based 
coordinates, but it should be noted that the resulting instabilities are generally moderate 
(not shown) and could be controlled by the use of slightly damping algorithms (e.g. SI 
decentering as in Laprise, 1995) or by the diffusive processes acting in a complete model. 

The stability of the SI scheme can be studied for non-isothermal linear atmospheric 
flows, using the "numerical analysis" method proposed by CBS83: The linearized equa- 
tions have first to be vertically discretized. For a given eigenmode of the horizontal 
operator, the actual and reference model operators can then be expressed as two matrices 
L and L* operating on the perturbation state-vector X. The equation of model evolution 
(fT2|) then writes: 

yt+At ^ ^ Y* (78) 

where the generalized state vector Y* is defined by (X*,X*~^*), and the amplification 

matrix A is given by: 

^ 2 At (Ix - At'L*y.\L - L*) (Ix - AtL*)"\lx + AtL*) ^ 



A 



(79) 



\ Ix Ox / 

where Ix and Ox are the identity and null operators in the state- vector space respectively. 

Noting r the largest modulus of A eigenvalues, the scheme is stable if F < 1 and unstable 
if F > 1. The eigenvector associated with F gives the vertical structure and polarisation 
of the most unstable mode. 

Our implementation of this "numerical analysis" method has been validated by com- 
parison with the results of the above analyses: the agreement was found very good (not 
shown). Here, the "numerical analysis" method is applied to quantify the impact of the 
change from (i to d in more realistic situations than the one assumed in the above anal- 
yses. Following SHB78, a realistic actual temperature profile is chosen: it consists in a 

tropospheric profile with a quasi-uniform dry static stability and an isothermal strato- 
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sphere above 200 hPa, as depicted in Fig. El (the surface temperature is set to 285 K). 
The vertical discretization of L and £* models follows BHBG95. The numerical analysis 
is performed for the same typical future NWP target values for (Ax, At) as in Fig. [H 
The stability is examined for values of T* varying in the interval [200 K, 320 K]. 

Figure [21 shows the growth rates obtained for 10, 20 and 30 regularly-spaced a levels, 
for the d and d variables. The theoretical disadvantage of d is confirmed: for high ver- 
tical resolutions, the scheme is unstable for almost every value of T*, with growth rates 
uncompatible with a practical use. The logarithm of the growth rate is found roughly pro- 
portional to the horizontal wavenumber /c, confirming that troubles linked to the presence 
of NL terms become more and more stringent when the horizontal resolution is increased. 
Moreover, reducing the timestep is found to be of no help for a given forecast range be- 
cause the logarithm of the growth rate is always roughly proportional to the time-step. 
Conversely, the use of variable d results in a stable scheme provided T* is chosen warm 
enough in the examined interval (and up to 440 K here, consistenly with analyses). 

Taking d as prognostic variable, the impact of the choice of the nonhydrostatic pressure 

variable ("P vs. P) can be examined numerically, still in the same context. It has been 

shown in section (HI why the system could be potentially unstable when vf^ deviates from 

TT*. To confirm and illustrate this statement, two values of tt* are chosen (813.25 hPa, and 

1213.25 hPa), while vf^ is assumed to be 1013.25 hPa, still for the same thermal profile 

and experimental context as above, and for 30 regularly-spaced a-levels. The solid lines in 

Figure ini show the growth rates obtained for the variable V (as explained in section (HI the 

growth rate has no dependency on vr* when the variable V is used, hence the corresponding 

growth rate for the variable V can be seen in Fig. O regardless of the value of vr*). The 

use of the V variable induces a global decrease in stability especially for /T*=1213.25 hPa 

and low values of T*. For high values of T*, an instability subsists for both values of vr*, 
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but it is weak (growth rates around 1.015), and the practical risk of using V for NWP 
is hence difficult to assess from this simplified theory and should be evaluated in a more 
realistic context. However, since there is no other potential advantage of using V instead 
of P, this question is probably not of primary interest. 

Finally, examination of the above algebra shows that the hybrid rj coordinate can 
also be responsible for instabilities when vf^ deviates from vr*, as it may e.g. be the case 
when the terrain is uniformly elevated. This is linked to the fact that in rj coordinate, all 
vertical operators {Q, 5, A/", d) deviate from their Sl-reference counterpart when Tf^ 7^ tt*, 
the resulting NL residuals potentially leading to instability. The discrepancy between the 
vertical metrics of C and £* models makes impossible the analysis for a general choice of 
A and B functions; hence, we use the "numerical analysis" method proposed by CBS83 
to examine the impact of using a hybrid coordinate instead of a pure terrain-following 
coordinate when lis 7^ tt*. 

The dashed lines in Figure 01 show the growth rates obtained for the variable V 
when a regularly flattening unstretched rj coordinate, defined by: Ti{rj) = —prj\n{rj)poo + 
?7 [1 + pln(?7)] TTg (with poo=1013.25 hPa, and p=0.25) is used with the same realistic ther- 
mal profile and discrepancies of ±200 hPa between vf^ and vr* as above. For high values of 
T*, the loss of stability is of similar magnitude as the one resulting from the choice of V 
as prognostic variable. It is noteworthy that the instability linked to the use of the hybrid 
coordinate in EE cannot be eliminated by choosing high values for T*, as it was the case 
for HPE system (Simmons and Temperton, 1997), hence there is indication that hybrid 
rj coordinate possesses a slight theoretical disadvantage from this point of view, but this 
may be not redhibitory in practice. In other respects, the rj coordinate has been claimed 
to bring some advantages compared to a coordinates, especially for the accuracy of the 

pressure-gradient force computation, and for the assimilation of high-level data (Simmons 
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and Burridge, 1981). Hence, according to the warning signal obtained here, the question 
of the relevance of the hybrid coordinate should legitimately have to be considered when 
implementing a real-case application in EE system with vr-type coordinate, but cannot be 
answered with a reasonable certainty at this stage. 

6 Comments 

The intrinsic sensivity of SI schemes stability to the choice of prognostic variables which 
has been demonstrated in this paper, is totally independent of the space discretization, 
since the analyses have been performed in the space-continuous context. As a conse- 
quence, any numerical model built with an intrinsically unstable variable will exhibit an 
unstable behaviour, regardless of the space-discretization and staggering choices that may 
be done. However, it has been shown that an intrisincally stable variable for height-based 
coordinates can become intrinsically unstable for mass-based coordinates because of the 
flow-dependent vertical metrics of this type of coordinates. Hence, it appears that for 
designing a SI numerical model in EE system, the choice of the prognostic variables can- 
not be made independently of the choice of the vertical coordinate, even at the level of 
space-continuous equations. 

One could wonder if this dependency is specific to the EE system or if it was already 
present in HPE systems, which are cast in vr-type coordinates. First it should be noted that 
HPE system does not offer a large latitude for the choice of prognostic variables: for the 
problems examined here, the prognostic variable for the vertically integrated continuity 
equation is the only relevant one: it can be set to be tt^ or g = ln(7rs). The above numerical 
method applied to the HPE system demonstrates the absence of sensitivity to the choice 

of these prognostic variables (not shown). For instance, the instabilities reported by 

23 



Simmons and Temperton (1997), which are the consequence of using a hybrid coordinate, 
develop almost identically for tt^ and q = ln(7rs) variables. The sensitivity of the stability 
to the choice of the set of prognostic variables discussed here is thus a specificity of fully 
elastic nonhydrostatic systems (anelastic systems have not been examined). 

7 Conclusion 

The stability of constant-coefficients SI schemes for the system of Euler Equations has 
been examined from an theoretical point of view, in deliberately simplified contexts to 
allow tractable analyses. The salient result of the study is that the stability can be 
dramatically affected by changes in the set of prognostic variables which are used to 
design the SI scheme, and this, independently of any space discretization. It appears that 
the choice of the two "nonhydrostatic" variables has to be carefully checked if an optimal 
stability is desired, especially for mass-based coordinates. Non-optimal choices generally 
result in growth rates incompatible with a practical NWP use, as experienced in BHBG95. 
This sensitivity of the SI scheme stability to the choice of prognostic variables was not 
present in the HPE system. 

If the SI reference surface-pressure deviates from its actual counterpart for mass-based 
coordinates, hybrid coordinates are found more unstable than the corresponding pure 
terrain-following coordinates, as reported previously for HPE system, but the practical 
consequences of this slight instability need to be evaluated in more realistic contexts. 

The analyses and results reported in this paper apply to constant-coefficients SI 

schemes, however, since only thermal and baric NL residuals are involved, the results 

extend identically to those SI schemes belonging to the class (ii) of the introduction for 

which the reference temperature and pressure are horizontally homogeneous (e.g. Thomas 
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et al., 1998; Qian et al, 1998). 

The practical interest of this theoretical study could be considered as quite limited, 
since only a very small part of the NL terms due to the discrepancy between the actual 
state X and the SI reference state X* has been considered. In this respect, the develop- 
ments presented here should rather be viewed as an endeavour to stress that a special care 
must be exerted in the choice of prognostic variables for building a constant-coefficients 
SI EE system, whatever type of coordinate is used. 

In a more general way, there is another practical interest to this study: as a predic- 
tive tool, it can serve for the validation of practical numerical applications: in a simi- 
lar way that analytically predicted stationnary orographic flows are commonly used for 
validating the relevance of space-discretisation schemes, the analytically predicted sta- 
bility can be used for a careful validation of time-discretisation schemes. Moreover, if a 
"numerical-analysis" tool is built, a complete set of cross-validations becomes possible be- 
tween analyses predictions, numerical-analyses diagnostics, and observed numerical-model 
behaviour. This may signiflcantly help to improve the detection of anomalous behaviour 
in a numerical model, or the prediction of the behaviour of alternative time discretizations 
(off-centered SI scheme, iterative treatments, etc.). 

Nevertheless, as stated above, this study unavoidably leads to an overestimation of 
the stability in comparison to the likely stability for real-case conditions because it retains 
only a very small part of all NL terms arising from the discrepancy {Ai — £*). Among the 
most important sources of NL residuals, the spatial variability of the orography has been 
neglected so far. Hence, there is no proof that the optimal set found here (e.g. q, V, d, 
and a coordinate) is still optimal, or even reasonnably stable, when a steep orography is 
introduced. Further examination of this topic would constitute an interesting extension 
for the present study. 
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Appendix A : List of Symbols 

vertical spatial operators in a coordinate: 
gX = £ {X/a')da' 
SX = {l/a)J^Xda' 
AfX = j^Xda 

XX = x 

dX = a{d/da)X 
miscellaneous symbols: 

V : Horizontal gradient operator along constant surfaces of the considered vertical 
coordinate. 

V : Horizontal wind vector in 3D framework. 

V : (d/dx) along constant levels of the considered vertical coordinate in the 2D verical 
plane domain. 

(d/dt) : Eulerian time-derivative 

(d/dt) : Lagrangian time-derivative 

D : horizontal wind divergence (V.m) 

g : gravitational acceleration 

p : true pressure 

vr : hydrostatic pressure 

Poo : absolute reference pressure (1013.25 hPa) 

R, Cp, : dry air thermodynamic constants 

u : horizontal wind component in the 2D framework 
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Appendix B : Stability analysis of the EE in regular Gal- 
Chen coordinate 



The formalism for the derivation of the constant-coefficients SI scheme analysed here 
basically follows Caya and Laprise (1999). Notations are standard ones and taken identical 
to this paper unless specified. The unstretched Gal-Chen coordinate is defined by: 

where Hd is the height of the domain, and ho is the height of the terrain. In the absence 
of orography {ho = 0) we thus have: 

C = z (81) 

The general framework is the same as in the above analyses, as depicted in section [21 

and the thermal discrepancy between X and X* states is still noted a = (T — T*)/T*. 

However, unlike for mass-based coordinates, the pressure variable q = ln(p/poo) (not 

to be confused with q = In^Hs) for mass-coordinates systems in the main part of the 

paper) needs to be defined in the whole space for both X and X* states. The hydrostatic 

equilibrium and stationnarity of these states implies: 

^dq ^ ^,d(f_ ^ _g_ 
dz dz R 

The SI time-discretized system is given by Eq. (46) - (50) of Caya and Laprise, 1999: 

^ + i?T*W* = -RT'Vq' (83) 

(84) 







dt dz 




oz 








dT' 
"dt ' 


RT* dq' 
Cp dt 


Op 


= -^T'{Vu + 


dw 
dz 


Cp 


'dq' 


RT* 


+ [ 


dw\ 

+ 


= 





(86) 

where x' = T — T* and q' = q — q*. This system is exact in the sense that no term has 

been neglected so far. As can be seen from the notation the left part is treated in a 
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centred implicit way while the right part is treated explicitly. For the analysis, a small 
pertubation (T, q) around the actual state (T, q) must be introduced and the complete 
system must be linearized around (T, q). Equating the total variables T and q yields: 

T = T'-aT* (87) 

The terms which are non-linear in (u, w, T, q) are then dropped in the above system 
(T' and q' cannot be considered as small however). The actual state X being a resting 
state, time derivatives are replaced by (d/dt). However, a special care must be taken to 
express (dq'/dt) in term of (dq/dt), because the vertical transport of (g — q*) is a linear 
term which must be retained. If an Eulerian scheme for the vertical advection is used, 
then: 

dq' dq a g 

where the advection term is treated explicitly. For a semi-Lagrangian vertical transport 
scheme, the assumptions of perfect solution for the displacement equation, and of perfect 
interpolators, consistent with the space-continuous context, lead after some manipulation 
based on Taylor series expansions in space, to the same result than for the Eulerian 
advection scheme as written in (fHn| . 

The linearization of the above system around the X state leads, after some algebraic 
manipulation to a form compatible with the analysis proposed in BOS. The C system 
then writes: 

^ = -RTVq (90) 

^ = If-RT^ (91) 
dt T dz ^ ^ 

g Cp ^ dw\ 



at Rr-ct[^"^al) ('^) 

dT RT dw\ , , 

V«-f-— . (93) 



dt \ dz 
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The SI reference system C* is shown to be formally identical, but with T replaced by T*. 
We have P = 4 and the normal modes of the time-continuous system are given by: 



where n = + 1/2)/ H and H = RT/g. It should be noted that the time-continuous 
normal modes of the C* system have a different vertical structure exp(n*2;) with n* = 
{iv + l/2)g/{RT*). This discrepancy between the height-scales for the vertical growth of 
C and C* time-continuous normal modes has some important consequences, as discussed 
in section O 

The analysis then proceeds as for mass-based coordinates by defining a numerical 
growth rate A. For a given geometry {k, u), the stability equation finally writes: 



where A_, A+, A, A, c, H, N are defined as in (IH^ - (f^ . 

The asymptotic growth rate is then given by the same equation (f75|) as for mass- 
based coordinates, hence the conditions for unconditional stability are the same as well. 
However, a major difference between Eqs. (f95|) and (f7i|l is the presence of the last complex 
term in the bracket of (|95|) . The solution of this equation for finite time-steps shows that 
this term, being complex, is responsible of an instability when a ^ 0. The instability is 
generally small, but is found to become quite significant when kH ^ z/ ^ 1, but, as stated 
above, it always disappears for long time-steps. 



Xj{x, z) = Xj exp{ikx) exp(?7.2;) j G (1, . . . 



(94) 




A;2(A+)A + nn(A+)A + nA=(A+-A) +eN\^A^A =0 (95) 
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List of Figures 

Fig. 1: Analytic growth rates (GR) for the simplified problem as a function of the 
nonlinearity parameter a. Thick line: asymptotic GR for variable d, solid line: asymptotic 
GR for variable d, dotted line: practical GR for variable d. dashed line: practical GR for 
variable d (see text for the meaning of practical GR). 

Fig. 2: Numerical growth rate (left axis) for the realistic thermal profile (dotted line 

- pressure on right axis), as a function of T*. From thin to thick lines: 10, 20, 30 levels, 
solid line: variable rf, dashed line: variable d. 

Fig. 3: Numerical growth rate (left axis) for the realistic thermal profile (dotted line 

- pressure on right function of T*. Thick lines: < = 1213.25 hPa, thin lines : 
vr* = 813.25 hPa, solid line: variable V with a coordinate, dashed line: variable V with rj 
coordinate. 
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Figure 1: Analytic growth rates (GR) for the simplified problem as a function of the 
nonlinearity parameter a. Thick line: asymptotic GR for variable d, solid line: asymptotic 
GR for variable d, dotted line: practical GR for variable d. dashed line: practical GR for 
variable d (see text for the meaning of practical GR). 
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Figure 2: Numerical growth rate (left axis) for the realistic thermal profile (dotted line 
- pressure on right function of T* . From thin to thick lines: 10, 20, 30 levels, 

solid line: variable d, dashed line: variable d. 
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Figure 3: Numerical growth rate (left axis) for the realistic thermal profile (dotted line 
- pressure on right axis), as a function of T*. Thick lines: tt* = 1213.25 hPa, thin lines : 
71* = 813.25 hPa, solid line: variable V with a coordinate, dashed line: variable V with r] 
coordinate. 
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